---
title: "Bats and human blood project: Sex-specific animal type analysis"
output:
  html_document:
    df_print: paged
  pdf_document: default
---

```{r setup, echo=FALSE}
knitr::opts_knit$set(root.dir = here::here())
# setwd(here::here())
```

# prepare & load packages
```{r}
# pathWD <- "/Volumes/5360/Biophysics/Experiments/dRT-DC/220301_Bob_DataAnalysis_BatsHumans/20230829_AnimalTypeTempSexAnalysis/"
pathWD <- "~/Documents/08_Pojekte/01_dRT-DC/220301_Bob_DataAnalysis_BatsHumans/20230829_AnimalTypeTempSexAnalysis"
setwd(pathWD)
# setwd(here::here())

pacman::p_load(tidyverse, lmerTest, car, emmeans, easystats, readr, ggh4x, ggplot2)

# load data ####
Bat1Table <- read_delim(file = "../Bat1Table_matProp.csv",
                        delim = ",",
                        col_names = TRUE)
Bat1Table$replicate <- paste0("NN", Bat1Table$index)

Bat2Table <- read_delim(file = "../Bat2Table_matProp.csv",
                        delim = ",",
                        col_names = TRUE)
Bat2Table$replicate <- paste0("NN", Bat2Table$index)

BatFLITable <- read_delim(file = "../BatFLITable_matProp.csv",
                          delim = ",",
                          col_names = TRUE)
BatFLITable$replicate <- paste("RA", BatFLITable$index, sep = "")

HumanTable <- read_delim(file = "../HumanTable_matProp.csv",
                         delim = ",",
                         col_names = TRUE)
HumanTable$replicate <- paste("H", HumanTable$index, sep = "")

# remove columns for consistency
Bat1Table <- select(Bat1Table, -bodyMass, -forearmLength, -tempAtCapture, -torporCapture, -torporBlood)

Bat2Table <- select(Bat2Table, -bodyMass, -forearmLength, -tempAtCapture, -torporCapture, -torporBlood)

# concatenate data tables
combTable <- rbind(Bat1Table,Bat2Table,BatFLITable,HumanTable)
remove(Bat1Table,Bat2Table,BatFLITable,HumanTable)

# combTable
# View(combTable)
```


# define factors & sample size tables
```{r}
# define factors ####
combTable$animalType <- as.factor(combTable$animalType)
combTable$replicate <- as.factor(combTable$replicate)
combTable$tempLevel <- as.factor(combTable$tempLevel)
combTable$sex <- as.factor(combTable$sex)
combTable$operator <- as.factor(combTable$operator)

# re-level species: "Human" "RousettusAegyptiacus" "NyctalusNoctula"
combTable$animalType <- relevel(combTable$animalType, ref = "RousettusAegyptiacus")
combTable$animalType <- relevel(combTable$animalType, ref = "Human")

# re-level temperature levels: "Warm" "RT" "Cold"
combTable$tempLevel <- relevel(combTable$tempLevel, ref = "RT")
combTable$tempLevel <- relevel(combTable$tempLevel, ref = "Warm")

# levels ####
levels(combTable$animalType)
levels(combTable$replicate)
levels(combTable$tempLevel)
levels(combTable$sex)
levels(combTable$operator)

# overview table ####
table(combTable$animalType)
table(combTable$tempLevel)
table(combTable$sex)

table(combTable$tempLevel, combTable$replicate)
table(combTable$animalType, combTable$replicate)
table(combTable$sex, combTable$replicate)
table(combTable$animalType, combTable$sex)
table(combTable$animalType, combTable$tempLevel)


```


# preparation
```{r}
levels.AnimalType <- levels(combTable$animalType)
levels.Temp <- levels(combTable$tempLevel)
levels.Sex <- levels(combTable$sex)

listOutcome <- c("area","tau1","tau2","deltaDefoTau1","deltaDefoTau2","eModulus","viscosity")
```

# density plots
```{r}
plots.Density <- purrr::map(listOutcome,
                           ~ ggplot(data = combTable,
                                    mapping = aes_(x = as.name(.x))) +
                             geom_density(size=1, alpha=0.3) +
                             ggtitle(paste(.x)))
plots.Density
```

# all models
```{r}
# formulas
listFormula.full <- purrr::map(listOutcome,
                               ~ as.formula(paste(.x, "~ animalType * tempLevel * sex + (1|replicate)"))
)

# calculate models
listModel.full <- purrr::map(listFormula.full,
                             ~ lmer(formula = .x,
                                    data = combTable,
                                    REML = TRUE,
                                    na.action = na.exclude)
)

# step function from lmerTest package
# remove interactions which are not significant
# but keep always the main effects
listModel.step <- purrr::map(listModel.full,
                             ~ lmerTest::step(.x,
                                              ddf = "Satterthwaite",
                                              alpha.fixed = 0.05,
                                              alpha.random = 0.1,
                                              reduce.fixed = TRUE,
                                              reduce.random = FALSE,
                                              keep = c("animalType", "tempLevel", "sex")
                             )
)

# extract final models
listModel.final <- purrr::map(listModel.step,
                              ~ get_model(.x)
)

```

# check all models
```{r, eval = FALSE}
# for (m in 1:length(listModel.full)) {
#   modelCheckOut <- check_model(listModel.full[[m]])
#   pdf(file = paste0("modelCheck_", listOutcome[m]),
#       width = 15,
#       height = 15)
#   print(modelCheckOut)
#   dev.off()
# }
```

# print outome
```{r}
for (m in 1:length(listModel.full)) {
  writeLines("------------------------------------------------------------")
  writeLines(paste("Outcome variable:", listOutcome[m]))
  # writeLines("------------------------------------------------------------")
  # writeLines(paste("Final model:", listModel.final[[m]]@call, "\n", "\n"))
  
  writeLines("-----Anova3-------------------------------------------------")
  # global p-values
  Anova(listModel.final[[m]], type="3") %>% print()
  
  writeLines("-----Summary------------------------------------------------")
  # summary (estimate) output
  summary(listModel.final[[m]], corr=FALSE) %>% print()
  
  writeLines("-----End----------------------------------------------------\n")
  
}
```

# fitted values
```{r}
listResults.fittedValues <- purrr::map(listModel.final,
                                       ~ data.frame(replicate = combTable$replicate,
                                                    animalType = combTable$animalType,
                                                    sex = combTable$sex,
                                                    tempLevel = combTable$tempLevel,
                                                    fitted = fitted.values(.x)))

# remove NA
listResults.fittedValues <- purrr::map(listResults.fittedValues,
                                       ~ na.omit(.x))

# remove all duplicates
listResults.fittedValues <- purrr::map(listResults.fittedValues,
                                       ~ .x[!duplicated(.x[c('replicate','animalType','tempLevel')]), ])
```


## mean over tempLevel & Warm condition
```{r}
library("dplyr")

listResults.fittedValues.mean <- purrr::map(listResults.fittedValues,
                                            ~ .x %>% group_by(replicate,animalType,sex) %>% summarise_at(vars(fitted), list(fitted.mean = mean)))

listResults.fittedValues.warm <- purrr::map(listResults.fittedValues,
                                            ~ subset(.x, tempLevel == "Warm"))
```


# plot preparation
```{r}
# library("ggh4x","tools","ggbeeswarm") # "tools" for capitalization
pacman::p_load("ggh4x","tools","ggbeeswarm") # "tools" for capitalization

# labels for animal type
# labels.AnimalType <- c("Human", expression(atop(italic(~Rousettus), italic(~aegyptiacus))), expression(atop(italic(~Nyctalus), italic(~noctula))))
# labels.AnimalType <- c("'Human'", "italic(~Rousettus)\n'aegyptiacus'", "'Nyctalus'\n'noctula'")
labels.AnimalType <- c("Homo\nsapiens", "Rousettus\naegyptiacus", "Nyctalus\nnoctula")
names(labels.AnimalType) <- levels.AnimalType

# labels.AnimalType.Short <- c("Human", expression(atop(italic(~Rou), italic(~aeg))), expression(atop(italic(~Nyc), italic(~noc))))
labels.AnimalType.Short <- c("Homo\nsap", "Rou\naeg", "Nyc\nnoc")
names(labels.AnimalType.Short) <- levels.AnimalType

labels.AnimalType.OneLine <- c(expression(italic("Homo sapiens")), expression(italic("Rousettus aegyptiacus")), expression(italic("Nyctalus noctula")))
names(labels.AnimalType.OneLine) <- levels.AnimalType

# labels for temp
labels.Temp <- c("37", "23", "10")
names(labels.Temp) <- levels.Temp
labels.TempUnit <- paste(labels.Temp, "°C")
names(labels.TempUnit) <- levels.Temp

# capital names for Sex
labels.Sex <- tools::toTitleCase(levels.Sex)
names(labels.Sex) <- levels.Sex

# axes titles
title.AnimalType <- "Species"
title.Temp <- "Temperature"
title.TempUnit <- paste(title.Temp, "(°C)")
title.Sex <- "Sex"

# labels for outcome variables
listLabelOutcome <- c(bquote("Cell size" ~ ("μm"^2)),
                      expression(tau["inlet"] ~ "(ms)"),
                      expression(tau["channel"] ~ "(ms)"),
                      expression(hat(d)["inlet"]),
                      expression(hat(d)["channel"]),
                      "Young's modulus (kPa)",
                      "Viscosity (Pa s)")

# color scheme
colorPalette.AnimalType <- c("#C21618","#F3BA29","#1F77BD")
colorPalette.AnimalTypeTemp <- rep(colorPalette.AnimalType, times = 3)
colorPalette.TempAnimalType <- rep(colorPalette.AnimalType, each = 3)
colorPalette.AnimalTypeSex <- c("#C21618","#E08A8B","#F3BA29","#F9DC94","#1F77BD","#8FBBDE")
colorPalette.SexAnimalType <- c("#C21618","#F3BA29","#1F77BD","#E08A8B","#F9DC94","#8FBBDE")

```


# animalType, results
```{r}
listResults.AnimalType <- purrr::map(listModel.final,
                                     ~ emmeans(.x,
                                               pairwise ~ animalType,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.AnimalType.means <- purrr::map(listResults.AnimalType,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.AnimalType.means <- purrr::map(listResults.AnimalType.means,
                                           ~ within(.x, {
                                             lower.SEM <- emmean - SE
                                             upper.SEM <- emmean + SE})
)

listResults.AnimalType.contrasts <- purrr::map(listResults.AnimalType,
                                               ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.AnimalType)
```


## plot preparation
```{r}
yMinFrac <- c(0.5, 0.2, 0.2, 0.5, 0.5, 0.7, 0.4)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.AnimalType.plotPara <- list()
n = nrow(listResults.AnimalType.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.AnimalType.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.AnimalType.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.AnimalType.contrasts[[m]] <- listResults.AnimalType.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1,1,2) + 0.05,
           xEnd     = c(2,3,3) - 0.05,
           yLineRel = (1 + c(0,1,0) ) / 2.8,
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots AnimalType
```{r fig.height=10/2.54, fig.width=9/2.54}
plots.AnimalType <- list()
tbl.plotResults.AnimalType <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.AnimalType.plotPara[[m]]$dataMax
  yLo <- listResults.AnimalType.plotPara[[m]]$yLo
  yHi <- listResults.AnimalType.plotPara[[m]]$yHi
  
  plots.AnimalType[[m]] <- ggplot(data = listResults.AnimalType.means[[m]])
  
  # define layout parameters
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    geom_bar(mapping = aes(x = animalType,
                           y = emmean,
                           fill = animalType),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    geom_errorbar(mapping = aes(x = animalType,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    scale_x_discrete(labels = labels.AnimalType)
  
  # add p-values
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    geom_segment(data = listResults.AnimalType.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.AnimalType.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  ymax <- listResults.AnimalType.means[[m]] %>% select(upper.CL) %>% max(., na.rm = TRUE)
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.AnimalType,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.AnimalType[[m]] <- plots.AnimalType[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalType", ".svg")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.AnimalType[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalType", ".png")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.AnimalType[[m]], scale = 1)
  
  # results table
  tbl.plotResults.AnimalType <- rbind(tbl.plotResults.AnimalType, data.frame(rep(listOutcome[m], times = nrow(listResults.AnimalType.means[[m]])), listResults.AnimalType.means[[m]]$animalType, listResults.AnimalType.means[[m]]$emmean, listResults.AnimalType.means[[m]]$SE, listResults.AnimalType.means[[m]]$upper.CL - listResults.AnimalType.means[[m]]$emmean))
}

colnames(tbl.plotResults.AnimalType) <- c("targetVar", "animalType", "mean", "SEM", "CI")
write.csv(tbl.plotResults.AnimalType, file.path(pathWD, paste0("tblPlotResults", "_AnimalType", ".csv")), row.names = FALSE)
plots.AnimalType

```


# animalType @ Warm, results
```{r}
listResults.AnimalType.warm <- purrr::map(listModel.final,
                                          ~ emmeans(.x,
                                                    pairwise ~ animalType | tempLevel,
                                                    lmer.df = "satterthwaite",
                                                    adjust = NULL,
                                                    lmerTest.limit = 1e6)
)

listResults.AnimalType.warm.means <- purrr::map(listResults.AnimalType.warm,
                                                ~ .x$emmeans %>% as.data.frame()
)
listResults.AnimalType.warm.means <- purrr::map(listResults.AnimalType.warm.means,
                                                ~ within(.x, {
                                                  lower.SEM <- emmean - SE
                                                  upper.SEM <- emmean + SE})
)

listResults.AnimalType.warm.contrasts <- purrr::map(listResults.AnimalType.warm,
                                                    ~ .x$contrasts %>% as.data.frame()
)

# gate on tempLevel == "Warm"
listResults.AnimalType.warm.means <- purrr::map(listResults.AnimalType.warm.means,
                                                ~ subset(.x, tempLevel == "Warm")
)

listResults.AnimalType.warm.contrasts <- purrr::map(listResults.AnimalType.warm.contrasts,
                                                    ~ subset(.x, tempLevel == "Warm")
)

rm(listResults.AnimalType.warm)
```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.AnimalType.warm.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.AnimalType.warm.contrasts[[m]] <- listResults.AnimalType.warm.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1,1,2) + 0.05,
           xEnd     = c(2,3,3) - 0.05,
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots
```{r fig.height=10/2.54, fig.width=9/2.54}
plots.AnimalType.warm.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.6, 0.5),
                                           max = c(46, 10, 10, 0.3, 0.3, 1.15, 3))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.AnimalType.warm.contrasts[[m]] <- listResults.AnimalType.warm.contrasts[[m]] %>%
    mutate(yStart = (0.85 + 0.08*c(0,1,0)) * (plots.AnimalType.warm.limits$max[m] - plots.AnimalType.warm.limits$min[m]) + plots.AnimalType.warm.limits$min[m],
           yEnd = yStart,
           yText = yStart + 0.015*(plots.AnimalType.warm.limits$max[m] - plots.AnimalType.warm.limits$min[m]))
}

plots.AnimalType.warm <- list()

for (m in 1:length(listModel.final)) {

  plots.AnimalType.warm[[m]] <- ggplot(data = combTable)

  # define layout parameters
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create point plot: replicate means
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    geom_beeswarm(data = listResults.fittedValues.warm[[m]],
                 mapping = aes(x = animalType,
                               y = fitted),
                 colour = "#C0C0C0",
                 size = 2,
                 cex = 1.6)
  
  # # create LMM plot: mean
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    geom_point(data = listResults.AnimalType.warm.means[[m]],
               mapping = aes(x = animalType,
                             y = emmean),
             stat = "identity", # important since we provide mean values already
             size = 3,
             colour = colorPalette.AnimalType)

  # # create LMM plot: error bars
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    geom_errorbar(data = listResults.AnimalType.warm.means[[m]],
                  mapping = aes(x = animalType,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = colorPalette.AnimalType,
                  position = position_dodge())
  
  # add p-values
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    geom_segment(data = listResults.AnimalType.warm.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.AnimalType.warm.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nested x-axis
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    scale_x_discrete(labels = labels.AnimalType)
  
  # add further layout parameter
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.AnimalType.warm.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.AnimalType,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.AnimalType.warm[[m]] <- plots.AnimalType.warm[[m]] +
    coord_cartesian(ylim = c(plots.AnimalType.warm.limits$min[m], plots.AnimalType.warm.limits$max[m]))
  
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalType", "_Warm", ".svg")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.AnimalType.warm[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalType", "_Warm", ".png")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.AnimalType.warm[[m]], scale = 1)
}

plots.AnimalType.warm

```


## plot prep: viscosity vs eModulus
```{r}
listResults.fittedValues.warm.viscEMod <- subset(listResults.fittedValues.warm[[which(listOutcome == "eModulus")]], select = c(replicate, animalType, sex, tempLevel))
listResults.fittedValues.warm.viscEMod <- listResults.fittedValues.warm.viscEMod %>%
  mutate(eMod = listResults.fittedValues.warm[[which(listOutcome == "eModulus")]]$fitted) %>%
  mutate(visc = listResults.fittedValues.warm[[which(listOutcome == "viscosity")]]$fitted)

listResults.AnimalType.warm.means.viscEMod <- subset(listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]], select = c(animalType, tempLevel))
listResults.AnimalType.warm.means.viscEMod <- listResults.AnimalType.warm.means.viscEMod %>%
  mutate(eMod.emmean = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$emmean,
         eMod.SE = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$SE,
         eMod.df = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$df,
         eMod.lower.CL = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$lower.CL,
         eMod.upper.CL = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$upper.CL,
         eMod.upper.SEM = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$upper.SEM,
         eMod.lower.SEM = listResults.AnimalType.warm.means[[which(listOutcome == "eModulus")]]$lower.SEM) %>%
  mutate(visc.emmean = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$emmean,
         visc.SE = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$SE,
         visc.df = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$df,
         visc.lower.CL = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$lower.CL,
         visc.upper.CL = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$upper.CL,
         visc.upper.SEM = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$upper.SEM,
         visc.lower.SEM = listResults.AnimalType.warm.means[[which(listOutcome == "viscosity")]]$lower.SEM)
```



## plots: viscosity vs eModulus
```{r fig.height=9.8/2.54, fig.width=9/2.54}
plots.AnimalType.warm.viscEMod.limits <- data.frame(min = c(0.5, 0.5),
                                                    max = c(0.9, 2.5))

plots.AnimalType.warm.viscEMod <- ggplot()

# define layout parameters
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.justification = c(0.02,.98), legend.position = c(0.02,.98))

# create point plot: replicate means
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  geom_point(data = listResults.fittedValues.warm.viscEMod,
             mapping = aes(x = eMod,
                           y = visc,
                           color = animalType),
             size = 2,
             alpha = 0.3)

# # create LMM plot: mean
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  geom_point(data = listResults.AnimalType.warm.means.viscEMod,
             mapping = aes(x = eMod.emmean,
                           y = visc.emmean,
                           color = animalType),
           stat = "identity", # important since we provide mean values already
           size = 3)

# # create LMM plot: error bars
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  geom_errorbar(data = listResults.AnimalType.warm.means.viscEMod,
                mapping = aes(x = eMod.emmean,
                              ymin = visc.lower.SEM,
                              ymax = visc.upper.SEM,
                              color = animalType),
                size = 0.75,
                width = 0.05*(plots.AnimalType.warm.viscEMod.limits$max[1] - plots.AnimalType.warm.viscEMod.limits$min[1]))
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  geom_errorbar(data = listResults.AnimalType.warm.means.viscEMod,
                mapping = aes(xmin = eMod.lower.SEM,
                              xmax = eMod.upper.SEM,
                              y = visc.emmean,
                              color = animalType),
                size = 0.75,
                width = 0.05*(plots.AnimalType.warm.viscEMod.limits$max[2] - plots.AnimalType.warm.viscEMod.limits$min[2]))

# add further layout parameter
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                     limits = c(0, plots.AnimalType.warm.limits$max[m]),
                     expand = expansion(mult = c(0.025, 0.025))) +
  labs(caption = paste("Models:", format(listFormula.full[[which(listOutcome == "eModulus")]]), "\n", format(listFormula.full[[which(listOutcome == "viscosity")]])),
       x = listLabelOutcome[which(listOutcome == "eModulus")],
       y = listLabelOutcome[which(listOutcome == "viscosity")])

# change legend
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  scale_colour_manual(name  = title.AnimalType,
                       breaks = names(labels.AnimalType.OneLine),
                       labels = labels.AnimalType.OneLine,
                       values = c("#E69F00", "#56B4E9", "#009E73")) 

# change coordinates = "zoom in"
plots.AnimalType.warm.viscEMod <- plots.AnimalType.warm.viscEMod +
  coord_cartesian(xlim = c(plots.AnimalType.warm.viscEMod.limits$min[1], plots.AnimalType.warm.viscEMod.limits$max[1]),
                  ylim = c(plots.AnimalType.warm.viscEMod.limits$min[2], plots.AnimalType.warm.viscEMod.limits$max[2]))


# save plots
ggsave(filename = file.path(pathWD, paste0("visc_vs_eMod", "_AnimalType", "_Warm", ".svg")),
       width = 9, height = 9.8, units = "cm", dpi = 600, device = "svg",
       plot = plots.AnimalType.warm.viscEMod, scale = 1)

ggsave(filename = file.path(pathWD, paste0("visc_vs_eMod", "_AnimalType", "_Warm", ".png")),
       width = 9, height = 9.8, units = "cm", dpi = 600, device = "png",
       plot = plots.AnimalType.warm.viscEMod, scale = 1)


plots.AnimalType.warm.viscEMod

```


# temp, results
```{r}
listResults.Temp <- purrr::map(listModel.final,
                               ~ emmeans(.x,
                                         pairwise ~ tempLevel,
                                         lmer.df = "satterthwaite",
                                         adjust = NULL,
                                         lmerTest.limit = 1e6)
)

listResults.Temp.means <- purrr::map(listResults.Temp,
                                     ~ .x$emmeans %>% as.data.frame()
)
listResults.Temp.means <- purrr::map(listResults.Temp.means,
                                     ~ within(.x, {
                                       lower.SEM <- emmean - SE
                                       upper.SEM <- emmean + SE})
)

listResults.Temp.contrasts <- purrr::map(listResults.Temp,
                                         ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.Temp)
```

## plot preparation
```{r}
yMinFrac <- c(0.8, 0.5, 0.5, 0.7, 0.7, 0.3, 0.1)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.Temp.plotPara <- list()
n = nrow(listResults.Temp.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.Temp.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.Temp.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.Temp.contrasts[[m]] <- listResults.Temp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1,1,2) + 0.05,
           xEnd     = c(2,3,3) - 0.05,
           yLineRel = (1 + c(0,1,0) ) / 2.8,
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots temp
```{r fig.height=10/2.54, fig.width=9/2.54}
plots.Temp <- list()
tbl.plotResults.Temp <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.Temp.plotPara[[m]]$dataMax
  yLo <- listResults.Temp.plotPara[[m]]$yLo
  yHi <- listResults.Temp.plotPara[[m]]$yHi
  
  plots.Temp[[m]] <- ggplot(data = listResults.Temp.means[[m]])
  
  # define layout parameters
  plots.Temp[[m]] <- plots.Temp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.Temp[[m]] <- plots.Temp[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.Temp[[m]] <- plots.Temp[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.Temp[[m]] <- plots.Temp[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.Temp[[m]] <- plots.Temp[[m]] +
    geom_segment(data = listResults.Temp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Temp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.Temp[[m]] <- plots.Temp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  ymax <- listResults.Temp.means[[m]] %>% select(upper.CL) %>% max(., na.rm = TRUE)
  plots.Temp[[m]] <- plots.Temp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.TempUnit,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.Temp[[m]] <- plots.Temp[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Temp", ".svg")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.Temp[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Temp", ".png")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.Temp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.Temp <- rbind(tbl.plotResults.Temp, data.frame(rep(listOutcome[m], times = nrow(listResults.Temp.means[[m]])), listResults.Temp.means[[m]]$tempLevel, listResults.Temp.means[[m]]$emmean, listResults.Temp.means[[m]]$SE, listResults.Temp.means[[m]]$upper.CL - listResults.Temp.means[[m]]$emmean))
}

colnames(tbl.plotResults.Temp) <- c("targetVar", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.Temp, file.path(pathWD, paste0("tblPlotResults", "_Temp", ".csv")), row.names = FALSE)
plots.Temp

```


# sex, results
```{r}
listResults.Sex <- purrr::map(listModel.final,
                               ~ emmeans(.x,
                                         pairwise ~ sex,
                                         lmer.df = "satterthwaite",
                                         adjust = NULL,
                                         lmerTest.limit = 1e6)
)

listResults.Sex.means <- purrr::map(listResults.Sex,
                                     ~ .x$emmeans %>% as.data.frame()
)
listResults.Sex.means <- purrr::map(listResults.Sex.means,
                                     ~ within(.x, {
                                       lower.SEM <- emmean - SE
                                       upper.SEM <- emmean + SE})
)

listResults.Sex.contrasts <- purrr::map(listResults.Sex,
                                         ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.Sex)
```

## plot preparation
```{r}
yMinFrac <- c(0.7, 0.5, 0.5, 0.6, 0.6, 0.8, 0.5)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.Sex.plotPara <- list()
n = nrow(listResults.Sex.contrasts[[1]])
for (m in 1:length(listModel.final)) {

  dataMax <- listResults.Sex.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.Sex.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo     = yLo,
                                               yHi     = yHi)
  
  listResults.Sex.contrasts[[m]] <- listResults.Sex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1) + 0.05,
           xEnd     = c(2) - 0.05,
           yLineRel = (1 + c(0) ) / 1.8,
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi)

}
```


## all plots sex
```{r fig.height=10/2.54, fig.width=9/2.54}
plots.Sex <- list()
tbl.plotResults.Sex <- data.frame(matrix(ncol = 4, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  dataMax <- listResults.Sex.plotPara[[m]]$dataMax
  yLo <- listResults.Sex.plotPara[[m]]$yLo
  yHi <- listResults.Sex.plotPara[[m]]$yHi
  
  plots.Sex[[m]] <- ggplot(data = listResults.Sex.means[[m]])
  
  # define layout parameters
  plots.Sex[[m]] <- plots.Sex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")
  
  # create plot
  plots.Sex[[m]] <- plots.Sex[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.Sex[[m]] <- plots.Sex[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.Sex[[m]] <- plots.Sex[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.Sex[[m]] <- plots.Sex[[m]] +
    geom_segment(data = listResults.Sex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Sex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.Sex[[m]] <- plots.Sex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  ymax <- listResults.Sex.means[[m]] %>% select(upper.CL) %>% max(., na.rm = TRUE)
  plots.Sex[[m]] <- plots.Sex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Sex,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.Sex[[m]] <- plots.Sex[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Sex", ".svg")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.Sex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Sex", ".png")),
         width = 9, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.Sex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.Sex <- rbind(tbl.plotResults.Sex, data.frame(rep(listOutcome[m], times = nrow(listResults.Sex.means[[m]])), listResults.Sex.means[[m]]$sex, listResults.Sex.means[[m]]$emmean, listResults.Sex.means[[m]]$SE, listResults.Sex.means[[m]]$upper.CL - listResults.Sex.means[[m]]$emmean))
}

colnames(tbl.plotResults.Sex) <- c("targetVar", "Sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.Sex, file.path(pathWD, paste0("tblPlotResults", "_Sex", ".csv")), row.names = FALSE)

plots.Sex

```


# animalType | temp, results
```{r}
listResults.AnimalTypeTemp <- purrr::map(listModel.final,
                                         ~ emmeans(.x,
                                                   pairwise ~ animalType | tempLevel,
                                                   lmer.df = "satterthwaite",
                                                   adjust = NULL,
                                                   lmerTest.limit = 1e6)
)

listResults.AnimalTypeTemp.means <- purrr::map(listResults.AnimalTypeTemp,
                                               ~ .x$emmeans %>% as.data.frame()
)
listResults.AnimalTypeTemp.means <- purrr::map(listResults.AnimalTypeTemp.means,
                                               ~ within(.x, {
                                                 lower.SEM <- emmean - SE
                                                 upper.SEM <- emmean + SE})
)

listResults.AnimalTypeTemp.contrasts <- purrr::map(listResults.AnimalTypeTemp,
                                                   ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.AnimalTypeTemp)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.AnimalTypeTemp.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.AnimalTypeTemp.contrasts[[m]] <- listResults.AnimalTypeTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Temp)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Temp)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots animalType | temp v2
```{r fig.height=10/2.54, fig.width=16/2.54}
plots.AnimalTypeTemp.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                          max = c(46, 10, 10, 0.3, 0.3, 1.9, 5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.AnimalTypeTemp.contrasts[[m]] <- listResults.AnimalTypeTemp.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.AnimalTypeTemp.limits$max[m] - plots.AnimalTypeTemp.limits$min[m]) + plots.AnimalTypeTemp.limits$min[m], times = length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.AnimalTypeTemp.limits$max[m] - plots.AnimalTypeTemp.limits$min[m]))
}

plots.AnimalTypeTemp <- list()
tbl.plotResults.AnimalTypeTemp <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.AnimalTypeTemp[[m]] <- ggplot(data = listResults.AnimalTypeTemp.means[[m]])
  
  # define layout parameters
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(cols = vars(tempLevel), labeller = labeller(tempLevel = labels.Temp))
  
  # create point plot: replicate means
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    geom_beeswarm(data = listResults.fittedValues[[m]],
                  mapping = aes(x = animalType,
                                y = fitted),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  
  # create LMM plot: mean
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    geom_point(mapping = aes(x = animalType,
                             y = emmean),
               stat = "identity", # important since we provide mean values already
               size = 3,
               colour = colorPalette.AnimalTypeTemp)
  
  # add LMM plot: error bars
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    geom_errorbar(mapping = aes(x = animalType,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = colorPalette.AnimalTypeTemp,
                  position = position_dodge())
  
  # add nested x-axis
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    scale_x_discrete(labels = labels.AnimalType.Short)
  
  # add p-values
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    geom_segment(data = listResults.AnimalTypeTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.AnimalTypeTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.AnimalTypeTemp.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.TempUnit,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.AnimalType,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.AnimalTypeTemp[[m]] <- plots.AnimalTypeTemp[[m]] +
    coord_cartesian(ylim = c(plots.AnimalTypeTemp.limits$min[m], plots.AnimalTypeTemp.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalTypeTemp", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.AnimalTypeTemp[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalTypeTemp", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.AnimalTypeTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.AnimalTypeTemp <- rbind(tbl.plotResults.AnimalTypeTemp, data.frame(rep(listOutcome[m], times = nrow(listResults.AnimalTypeTemp.means[[1]])), listResults.AnimalTypeTemp.means[[m]]$animalType, listResults.AnimalTypeTemp.means[[m]]$tempLevel,  listResults.AnimalTypeTemp.means[[m]]$emmean,  listResults.AnimalTypeTemp.means[[m]]$SE, listResults.AnimalTypeTemp.means[[m]]$upper.CL - listResults.AnimalTypeTemp.means[[m]]$emmean))
}

colnames(tbl.plotResults.AnimalTypeTemp) <- c("targetVar", "animalType", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.AnimalTypeTemp, file.path(pathWD, paste0("tblPlotResults", "_AnimalTypeTemp", ".csv")), row.names = FALSE)
plots.AnimalTypeTemp
```


## plot prep: viscosity vs eModulus
```{r}
listResults.fittedValues.viscEMod <- subset(listResults.fittedValues[[which(listOutcome == "eModulus")]], select = c(replicate, animalType, sex, tempLevel))
listResults.fittedValues.viscEMod <- listResults.fittedValues.viscEMod %>%
  mutate(eMod = listResults.fittedValues[[which(listOutcome == "eModulus")]]$fitted) %>%
  mutate(visc = listResults.fittedValues[[which(listOutcome == "viscosity")]]$fitted)

listResults.AnimalTypeTemp.means.viscEMod <- subset(listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]], select = c(animalType, tempLevel))
listResults.AnimalTypeTemp.means.viscEMod <- listResults.AnimalTypeTemp.means.viscEMod %>%
  mutate(eMod.emmean = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$emmean,
         eMod.SE = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$SE,
         eMod.df = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$df,
         eMod.lower.CL = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$lower.CL,
         eMod.upper.CL = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$upper.CL,
         eMod.upper.SEM = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$upper.SEM,
         eMod.lower.SEM = listResults.AnimalTypeTemp.means[[which(listOutcome == "eModulus")]]$lower.SEM) %>%
  mutate(visc.emmean = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$emmean,
         visc.SE = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$SE,
         visc.df = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$df,
         visc.lower.CL = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$lower.CL,
         visc.upper.CL = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$upper.CL,
         visc.upper.SEM = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$upper.SEM,
         visc.lower.SEM = listResults.AnimalTypeTemp.means[[which(listOutcome == "viscosity")]]$lower.SEM)
```



## plots: viscosity vs eModulus
```{r fig.height=9.8/2.54, fig.width=18/2.54}
plots.AnimalTypeTemp.viscEMod.limits <- data.frame(min = c(0.5, 0.5),
                                                    max = c(1.62, 4.4))

plots.AnimalTypeTemp.viscEMod <- ggplot(data = listResults.AnimalTypeTemp.means.viscEMod)

# define layout parameters
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.justification = c(0.02,.98), legend.position = c(0.02,.98)) +
    facet_grid(cols = vars(tempLevel), labeller = labeller(tempLevel = labels.Temp))

# create point plot: replicate means
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  geom_point(data = listResults.fittedValues.viscEMod,
             mapping = aes(x = eMod,
                           y = visc,
                           color = animalType),
             size = 2,
             alpha = 0.2)

# # create LMM plot: mean
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  geom_point(data = listResults.AnimalTypeTemp.means.viscEMod,
             mapping = aes(x = eMod.emmean,
                           y = visc.emmean,
                           color = animalType),
           stat = "identity", # important since we provide mean values already
           size = 1.5)

# # create LMM plot: error bars
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  geom_errorbar(data = listResults.AnimalTypeTemp.means.viscEMod,
                mapping = aes(x = eMod.emmean,
                              ymin = visc.lower.SEM,
                              ymax = visc.upper.SEM,
                              color = animalType),
                size = 0.5,
                width = 0.03*(plots.AnimalTypeTemp.viscEMod.limits$max[1] - plots.AnimalTypeTemp.viscEMod.limits$min[1]))
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  geom_errorbar(data = listResults.AnimalTypeTemp.means.viscEMod,
                mapping = aes(xmin = eMod.lower.SEM,
                              xmax = eMod.upper.SEM,
                              y = visc.emmean,
                              color = animalType),
                size = 0.5,
                width = 0.03*(plots.AnimalTypeTemp.viscEMod.limits$max[2] - plots.AnimalTypeTemp.viscEMod.limits$min[2]))

# add further layout parameter
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                     limits = c(0, plots.AnimalTypeTemp.limits$max[m]),
                     expand = expansion(mult = c(0.025, 0.025))) +
  labs(subtitle = title.Temp,
       caption = paste("Models:", format(listFormula.full[[which(listOutcome == "eModulus")]]), "\n", format(listFormula.full[[which(listOutcome == "viscosity")]])),
       x = listLabelOutcome[which(listOutcome == "eModulus")],
       y = listLabelOutcome[which(listOutcome == "viscosity")])

# change legend
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  scale_colour_manual(name  = title.AnimalType,
                       breaks = names(labels.AnimalType.OneLine),
                       labels = labels.AnimalType.OneLine,
                       values = c("#E69F00", "#56B4E9", "#009E73")) 

# change coordinates = "zoom in"
plots.AnimalTypeTemp.viscEMod <- plots.AnimalTypeTemp.viscEMod +
  coord_cartesian(xlim = c(plots.AnimalTypeTemp.viscEMod.limits$min[1], plots.AnimalTypeTemp.viscEMod.limits$max[1]),
                  ylim = c(plots.AnimalTypeTemp.viscEMod.limits$min[2], plots.AnimalTypeTemp.viscEMod.limits$max[2]))


# save plots
ggsave(filename = file.path(pathWD, paste0("visc_vs_eMod", "_AnimalTypeTemp", ".svg")),
       width = 18, height = 9.8, units = "cm", dpi = 600, device = "svg",
       plot = plots.AnimalTypeTemp.viscEMod, scale = 1)

ggsave(filename = file.path(pathWD, paste0("visc_vs_eMod", "_AnimalTypeTemp", ".png")),
       width = 18, height = 9.8, units = "cm", dpi = 600, device = "png",
       plot = plots.AnimalTypeTemp.viscEMod, scale = 1)


plots.AnimalTypeTemp.viscEMod

```


# temp | animalType, results
```{r}
listResults.TempAnimalType <- purrr::map(listModel.final,
                                         ~ emmeans(.x,
                                                   pairwise ~ tempLevel | animalType,
                                                   lmer.df = "satterthwaite",
                                                   adjust = NULL,
                                                   lmerTest.limit = 1e6)
)

listResults.TempAnimalType.means <- purrr::map(listResults.TempAnimalType,
                                               ~ .x$emmeans %>% as.data.frame()
)
listResults.TempAnimalType.means <- purrr::map(listResults.TempAnimalType.means,
                                               ~ within(.x, {
                                                 lower.SEM <- emmean - SE
                                                 upper.SEM <- emmean + SE})
)

listResults.TempAnimalType.contrasts <- purrr::map(listResults.TempAnimalType,
                                                   ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.TempAnimalType)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.TempAnimalType.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.TempAnimalType.contrasts[[m]] <- listResults.TempAnimalType.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.AnimalType)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.AnimalType)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots temp | animalType v2
```{r fig.height=10/2.54, fig.width=16/2.54}
plots.TempAnimalType.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                          max = c(46, 10, 10, 0.3, 0.3, 1.9, 5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.TempAnimalType.contrasts[[m]] <- listResults.TempAnimalType.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.TempAnimalType.limits$max[m] - plots.TempAnimalType.limits$min[m]) + plots.TempAnimalType.limits$min[m], times = length(levels.AnimalType)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.TempAnimalType.limits$max[m] - plots.TempAnimalType.limits$min[m]))
}

plots.TempAnimalType <- list()
tbl.plotResults.TempAnimalType <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.TempAnimalType[[m]] <- ggplot(data = listResults.TempAnimalType.means[[m]])
  
  # define layout parameters
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(cols = vars(animalType), labeller = labeller(animalType = labels.AnimalType))
  
  # create point plot: replicate means
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    geom_beeswarm(data = listResults.fittedValues[[m]],
                  mapping = aes(x = tempLevel,
                                y = fitted),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  
  # create LMM plot: mean
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    geom_point(mapping = aes(x = tempLevel,
                             y = emmean),
               stat = "identity", # important since we provide mean values already
               size = 3,
               colour = colorPalette.TempAnimalType)
  
  # add LMM plot: error bars
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = colorPalette.TempAnimalType,
                  position = position_dodge())
  
  # add nested x-axis
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    geom_segment(data = listResults.TempAnimalType.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempAnimalType.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.TempAnimalType.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.AnimalType,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.TempUnit,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.TempAnimalType[[m]] <- plots.TempAnimalType[[m]] +
    coord_cartesian(ylim = c(plots.TempAnimalType.limits$min[m], plots.TempAnimalType.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempAnimalType", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempAnimalType[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempAnimalType", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.TempAnimalType[[m]], scale = 1)
  
  # results table
  tbl.plotResults.TempAnimalType <- rbind(tbl.plotResults.TempAnimalType, data.frame(rep(listOutcome[m], times = nrow(listResults.TempAnimalType.means[[1]])), listResults.TempAnimalType.means[[m]]$tempLevel, listResults.TempAnimalType.means[[m]]$animalType, listResults.TempAnimalType.means[[m]]$emmean, listResults.TempAnimalType.means[[m]]$SE, listResults.TempAnimalType.means[[m]]$upper.CL - listResults.TempAnimalType.means[[m]]$emmean))
}

colnames(tbl.plotResults.TempAnimalType) <- c("targetVar", "tempLevel", "animalType", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TempAnimalType, file.path(pathWD, paste0("tblPlotResults", "_TempAnimalType", ".csv")), row.names = FALSE)
plots.TempAnimalType

```


# sex | animal type, results
```{r}
listResults.SexAnimalType <- purrr::map(listModel.final,
                               ~ emmeans(.x,
                                         pairwise ~ sex | animalType,
                                         lmer.df = "satterthwaite",
                                         adjust = NULL,
                                         lmerTest.limit = 1e6)
)

listResults.SexAnimalType.means <- purrr::map(listResults.SexAnimalType,
                                     ~ .x$emmeans %>% as.data.frame()
)
listResults.SexAnimalType.means <- purrr::map(listResults.SexAnimalType.means,
                                               ~ within(.x, {
                                                 lower.SEM <- emmean - SE
                                                 upper.SEM <- emmean + SE})
)

listResults.SexAnimalType.contrasts <- purrr::map(listResults.SexAnimalType,
                                         ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.SexAnimalType)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.SexAnimalType.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  
  listResults.SexAnimalType.contrasts[[m]] <- listResults.SexAnimalType.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1) + 0.05, times = length(levels.AnimalType)),
           xEnd     = rep(c(2) - 0.05, times = length(levels.AnimalType)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots sex | animal type
```{r fig.height=10/2.54, fig.width=14/2.54}
plots.SexAnimalType.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.0, 0.0),
                                             max = c(45, 10, 10, 0.3, 0.3, 1.75, 4.3))

# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.SexAnimalType.contrasts[[m]] <- listResults.SexAnimalType.contrasts[[m]] %>%
    mutate(yStart = rep((0.92 + 0.08*c(0)) * (plots.SexAnimalType.limits$max[m] - plots.SexAnimalType.limits$min[m]) + plots.SexAnimalType.limits$min[m], times = length(levels.AnimalType)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.SexAnimalType.limits$max[m] - plots.SexAnimalType.limits$min[m]))
}

plots.SexAnimalType <- list()
tbl.plotResults.SexAnimalType <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.SexAnimalType[[m]] <- ggplot(data = listResults.SexAnimalType.means[[m]])
  
  # define layout parameters
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 12, angle = 90),
          strip.text.x = element_text(size = 12),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(animalType), labeller = labeller(animalType = labels.AnimalType))
  
  # create plot
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             fill = colorPalette.AnimalTypeSex,
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    geom_segment(data = listResults.SexAnimalType.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.SexAnimalType.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  # plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
  #   scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.SexAnimalType.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.AnimalType,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Sex,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.SexAnimalType[[m]] <- plots.SexAnimalType[[m]] +
    coord_cartesian(ylim = c(plots.SexAnimalType.limits$min[m], plots.SexAnimalType.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexAnimalType", ".svg")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.SexAnimalType[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexAnimalType", ".png")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.SexAnimalType[[m]], scale = 1)
  
  # results table
  tbl.plotResults.SexAnimalType <- rbind(tbl.plotResults.SexAnimalType, data.frame(rep(listOutcome[m], times = nrow(listResults.SexAnimalType.means[[1]])), listResults.SexAnimalType.means[[m]]$sex,  listResults.SexAnimalType.means[[m]]$animalType, listResults.SexAnimalType.means[[m]]$emmean, listResults.SexAnimalType.means[[m]]$SE, listResults.SexAnimalType.means[[m]]$upper.CL - listResults.SexAnimalType.means[[m]]$emmean))
}

colnames(tbl.plotResults.SexAnimalType) <- c("targetVar", "sex", "animalType", "mean", "SEM", "CI")
write.csv(tbl.plotResults.SexAnimalType, file.path(pathWD, paste0("tblPlotResults", "_SexAnimalType", ".csv")), row.names = FALSE)

plots.SexAnimalType

```


# animal type | sex, results
```{r}
listResults.AnimalTypeSex <- purrr::map(listModel.final,
                               ~ emmeans(.x,
                                         pairwise ~ animalType | sex,
                                         lmer.df = "satterthwaite",
                                         adjust = NULL,
                                         lmerTest.limit = 1e6)
)

listResults.AnimalTypeSex.means <- purrr::map(listResults.AnimalTypeSex,
                                     ~ .x$emmeans %>% as.data.frame()
)
listResults.AnimalTypeSex.means <- purrr::map(listResults.AnimalTypeSex.means,
                                               ~ within(.x, {
                                                 lower.SEM <- emmean - SE
                                                 upper.SEM <- emmean + SE})
)

listResults.AnimalTypeSex.contrasts <- purrr::map(listResults.AnimalTypeSex,
                                         ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.AnimalTypeSex)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.AnimalTypeSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.AnimalTypeSex.contrasts[[m]] <- listResults.AnimalTypeSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Sex)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Sex)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots animal type | sex v2
```{r fig.height=10/2.54, fig.width=14/2.54}
plots.AnimalTypeSex.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                         max = c(45, 10, 10, 0.3, 0.3, 2, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.AnimalTypeSex.contrasts[[m]] <- listResults.AnimalTypeSex.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.AnimalTypeSex.limits$max[m] - plots.AnimalTypeSex.limits$min[m]) + plots.AnimalTypeSex.limits$min[m], times = length(levels.Sex)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.AnimalTypeSex.limits$max[m] - plots.AnimalTypeSex.limits$min[m]))
}

plots.AnimalTypeSex <- list()
tbl.plotResults.AnimalTypeSex <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.AnimalTypeSex[[m]] <- ggplot(data = listResults.AnimalTypeSex.means[[m]])
  
  # define layout parameters
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(sex), labeller = labeller(sex = labels.Sex))
  
  # create point plot: replicate means
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    geom_beeswarm(data = listResults.fittedValues[[m]],
                  mapping = aes(x = animalType,
                                y = fitted),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  
  # create LMM plot: mean
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    geom_point(mapping = aes(x = animalType,
                             y = emmean),
               stat = "identity", # important since we provide mean values already
               size = 3,
               colour = "#1f77bd")
  
  # add LMM plot: error bars
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    geom_errorbar(mapping = aes(x = animalType,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = "#1f77bd",
                  position = position_dodge())
  
  # add nested x-axis
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    scale_x_discrete(labels = labels.AnimalType)
  
  # add p-values
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    geom_segment(data = listResults.AnimalTypeSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.AnimalTypeSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.AnimalTypeSex.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Sex,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.AnimalType,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.AnimalTypeSex[[m]] <- plots.AnimalTypeSex[[m]] +
    coord_cartesian(ylim = c(plots.AnimalTypeSex.limits$min[m], plots.AnimalTypeSex.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalTypeSex", ".svg")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.AnimalTypeSex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalTypeSex", ".png")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.AnimalTypeSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.AnimalTypeSex <- rbind(tbl.plotResults.AnimalTypeSex, data.frame(rep(listOutcome[m], times = nrow(listResults.AnimalTypeSex.means[[1]])), listResults.AnimalTypeSex.means[[m]]$animalType, listResults.AnimalTypeSex.means[[m]]$sex,  listResults.AnimalTypeSex.means[[m]]$emmean,  listResults.AnimalTypeSex.means[[m]]$SE, listResults.AnimalTypeSex.means[[m]]$upper.CL - listResults.AnimalTypeSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.AnimalTypeSex) <- c("targetVar", "animalType", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.AnimalTypeSex, file.path(pathWD, paste0("tblPlotResults", "_AnimalTypeSex", ".csv")), row.names = FALSE)

plots.AnimalTypeSex

```


# temp | sex, results
```{r}
listResults.TempSex <- purrr::map(listModel.final,
                               ~ emmeans(.x,
                                         pairwise ~ tempLevel | sex,
                                         lmer.df = "satterthwaite",
                                         adjust = NULL,
                                         lmerTest.limit = 1e6)
)

listResults.TempSex.means <- purrr::map(listResults.TempSex,
                                     ~ .x$emmeans %>% as.data.frame()
)
listResults.TempSex.means <- purrr::map(listResults.TempSex.means,
                                         ~ within(.x, {
                                           lower.SEM <- emmean - SE
                                           upper.SEM <- emmean + SE})
)

listResults.TempSex.contrasts <- purrr::map(listResults.TempSex,
                                         ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.TempSex)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.TempSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.TempSex.contrasts[[m]] <- listResults.TempSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Sex)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Sex)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots temp | sex
```{r fig.height=10/2.54, fig.width=14/2.54}
plots.TempSex.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                   max = c(45, 10, 10, 0.3, 0.3, 2, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.TempSex.contrasts[[m]] <- listResults.TempSex.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.TempSex.limits$max[m] - plots.TempSex.limits$min[m]) + plots.TempSex.limits$min[m], times = length(levels.Sex)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.TempSex.limits$max[m] - plots.TempSex.limits$min[m]))
}

plots.TempSex <- list()
tbl.plotResults.TempSex <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.TempSex[[m]] <- ggplot(data = listResults.TempSex.means[[m]])
  
  # define layout parameters
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(sex), labeller = labeller(sex = labels.Sex))
  
  # create plot
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    geom_segment(data = listResults.TempSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.TempSex.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Sex,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.TempUnit,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.TempSex[[m]] <- plots.TempSex[[m]] +
    coord_cartesian(ylim = c(plots.TempSex.limits$min[m], plots.TempSex.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempSex", ".svg")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempSex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempSex", ".png")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.TempSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.TempSex <- rbind(tbl.plotResults.TempSex, data.frame(rep(listOutcome[m], times = nrow(listResults.TempSex.means[[1]])), listResults.TempSex.means[[m]]$tempLevel, listResults.TempSex.means[[m]]$sex,  listResults.TempSex.means[[m]]$emmean,  listResults.TempSex.means[[m]]$SE, listResults.TempSex.means[[m]]$upper.CL - listResults.TempSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.TempSex) <- c("targetVar", "tempLevel", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TempSex, file.path(pathWD, paste0("tblPlotResults", "_TempSex", ".csv")), row.names = FALSE)

plots.TempSex

```


# sex | temp, results
```{r}
listResults.SexTemp <- purrr::map(listModel.final,
                                         ~ emmeans(.x,
                                                   pairwise ~ sex | tempLevel,
                                                   lmer.df = "satterthwaite",
                                                   adjust = NULL,
                                                   lmerTest.limit = 1e6)
)

listResults.SexTemp.means <- purrr::map(listResults.SexTemp,
                                               ~ .x$emmeans %>% as.data.frame()
)
listResults.SexTemp.means <- purrr::map(listResults.SexTemp.means,
                                         ~ within(.x, {
                                           lower.SEM <- emmean - SE
                                           upper.SEM <- emmean + SE})
)

listResults.SexTemp.contrasts <- purrr::map(listResults.SexTemp,
                                                   ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.SexTemp)
```

## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.SexTemp.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.SexTemp.contrasts[[m]] <- listResults.SexTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1) + 0.05, times = length(levels.Temp)),
           xEnd     = rep(c(2) - 0.05, times = length(levels.Temp)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots sex | temp
```{r fig.height=10/2.54, fig.width=14/2.54}
plots.SexTemp.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                          max = c(45, 10, 10, 0.3, 0.3, 2, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.SexTemp.contrasts[[m]] <- listResults.SexTemp.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(1)) * (plots.SexTemp.limits$max[m] - plots.SexTemp.limits$min[m]) + plots.SexTemp.limits$min[m], times = length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.SexTemp.limits$max[m] - plots.SexTemp.limits$min[m]))
}

plots.SexTemp <- list()
tbl.plotResults.SexTemp <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.SexTemp[[m]] <- ggplot(data = listResults.SexTemp.means[[m]])
  
  # define layout parameters
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(cols = vars(tempLevel), labeller = labeller(tempLevel = labels.Temp))
  
  # create plot
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    geom_segment(data = listResults.SexTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.SexTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.SexTemp.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.TempUnit,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.Sex,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.SexTemp[[m]] <- plots.SexTemp[[m]] +
    coord_cartesian(ylim = c(plots.SexTemp.limits$min[m], plots.SexTemp.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexTemp", ".svg")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.SexTemp[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexTemp", ".png")),
         width = 14, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.SexTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.SexTemp <- rbind(tbl.plotResults.SexTemp, data.frame(rep(listOutcome[m], times = nrow(listResults.SexTemp.means[[1]])), listResults.SexTemp.means[[m]]$sex, listResults.SexTemp.means[[m]]$tempLevel,  listResults.SexTemp.means[[m]]$emmean,  listResults.SexTemp.means[[m]]$SE, listResults.SexTemp.means[[m]]$upper.CL - listResults.SexTemp.means[[m]]$emmean))
}

colnames(tbl.plotResults.SexTemp) <- c("targetVar", "sex", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.SexTemp, file.path(pathWD, paste0("tblPlotResults", "_SexTemp", ".csv")), row.names = FALSE)

plots.SexTemp

```


# animal type | temp * sex, results
```{r}
listResults.AnimalTypeTempSex <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ animalType | tempLevel * sex,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.AnimalTypeTempSex.means <- purrr::map(listResults.AnimalTypeTempSex,
                                                   ~ .x$emmeans %>% as.data.frame()
)
listResults.AnimalTypeTempSex.means <- purrr::map(listResults.AnimalTypeTempSex.means,
                                                   ~ within(.x, {
                                                     lower.SEM <- emmean - SE
                                                     upper.SEM <- emmean + SE})
)

listResults.AnimalTypeTempSex.contrasts <- purrr::map(listResults.AnimalTypeTempSex,
                                                       ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.AnimalTypeTempSex)
```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.AnimalTypeTempSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.AnimalTypeTempSex.contrasts[[m]] <- listResults.AnimalTypeTempSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Temp) * length(levels.Sex)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Temp) * length(levels.Sex)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots animal type | temp * sex
```{r fig.height=16/2.54, fig.width=16/2.54}
plots.AnimalTypeTempSex.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                          max = c(45, 10, 10, 0.3, 0.3, 2, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.AnimalTypeTempSex.contrasts[[m]] <- listResults.AnimalTypeTempSex.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.AnimalTypeTempSex.limits$max[m] - plots.AnimalTypeTempSex.limits$min[m]) + plots.AnimalTypeTempSex.limits$min[m], times = length(levels.Temp) * length(levels.Sex)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.AnimalTypeTempSex.limits$max[m] - plots.AnimalTypeTempSex.limits$min[m]))
}

plots.AnimalTypeTempSex <- list()
tbl.plotResults.AnimalTypeTempSex <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.AnimalTypeTempSex[[m]] <- ggplot(data = listResults.AnimalTypeTempSex.means[[m]])
  
  # define layout parameters
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(rows = vars(sex), cols = vars(tempLevel), labeller = labeller(sex = labels.Sex, tempLevel = labels.Temp))
  
  # create plot
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    geom_bar(mapping = aes(x = animalType,
                           y = emmean,
                           fill = animalType),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    geom_errorbar(mapping = aes(x = animalType,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    scale_x_discrete(labels = labels.AnimalType.Short)
  
  # add p-values
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    geom_segment(data = listResults.AnimalTypeTempSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.AnimalTypeTempSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.AnimalTypeTempSex.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.TempUnit,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.AnimalType,
         y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.AnimalTypeTempSex[[m]] <- plots.AnimalTypeTempSex[[m]] +
    coord_cartesian(ylim = c(plots.AnimalTypeTempSex.limits$min[m], plots.AnimalTypeTempSex.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalTypeTempSex", ".svg")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.AnimalTypeTempSex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_AnimalTypeTempSex", ".png")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.AnimalTypeTempSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.AnimalTypeTempSex <- rbind(tbl.plotResults.AnimalTypeTempSex, data.frame(rep(listOutcome[m], times = nrow(listResults.AnimalTypeTempSex.means[[1]])), listResults.AnimalTypeTempSex.means[[m]]$animalType, listResults.AnimalTypeTempSex.means[[m]]$tempLevel,  listResults.AnimalTypeTempSex.means[[m]]$sex, listResults.AnimalTypeTempSex.means[[m]]$emmean, listResults.AnimalTypeTempSex.means[[m]]$SE, listResults.AnimalTypeTempSex.means[[m]]$upper.CL - listResults.AnimalTypeTempSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.AnimalTypeTempSex) <- c("targetVar", "animalType", "tempLevel", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.AnimalTypeTempSex, file.path(pathWD, paste0("tblPlotResults", "_AnimalTypeTempSex", ".csv")), row.names = FALSE)

plots.AnimalTypeTempSex

```


# temp | animal type * sex, results
```{r}
listResults.TempAnimalTypeSex <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ tempLevel | animalType * sex,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.TempAnimalTypeSex.means <- purrr::map(listResults.TempAnimalTypeSex,
                                             ~ .x$emmeans %>% as.data.frame()
)
listResults.TempAnimalTypeSex.means <- purrr::map(listResults.TempAnimalTypeSex.means,
                                                   ~ within(.x, {
                                                     lower.SEM <- emmean - SE
                                                     upper.SEM <- emmean + SE})
)

listResults.TempAnimalTypeSex.contrasts <- purrr::map(listResults.TempAnimalTypeSex,
                                                 ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.TempAnimalTypeSex)
```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.TempAnimalTypeSex.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.TempAnimalTypeSex.contrasts[[m]] <- listResults.TempAnimalTypeSex.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.AnimalType) * length(levels.Sex)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.AnimalType) * length(levels.Sex)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots temp | animal type * sex
```{r fig.height=16/2.54, fig.width=16/2.54}
plots.TempAnimalTypeSex.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                          max = c(45, 10, 10, 0.3, 0.3, 2, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.TempAnimalTypeSex.contrasts[[m]] <- listResults.TempAnimalTypeSex.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.TempAnimalTypeSex.limits$max[m] - plots.TempAnimalTypeSex.limits$min[m]) + plots.TempAnimalTypeSex.limits$min[m], times = length(levels.AnimalType) * length(levels.Sex)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.TempAnimalTypeSex.limits$max[m] - plots.TempAnimalTypeSex.limits$min[m]))
}

plots.TempAnimalTypeSex <- list()
tbl.plotResults.TempAnimalTypeSex <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.TempAnimalTypeSex[[m]] <- ggplot(data = listResults.TempAnimalTypeSex.means[[m]])
  
  # define layout parameters
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(rows = vars(sex), cols = vars(animalType), labeller = labeller(sex = labels.Sex, animalType = labels.AnimalType))
  
  # create plot
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    geom_segment(data = listResults.TempAnimalTypeSex.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempAnimalTypeSex.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  ymax <- listResults.TempAnimalTypeSex.means[[m]] %>% select(upper.CL) %>% max(., na.rm = TRUE)
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.TempAnimalTypeSex.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.AnimalType,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = paste(title.TempUnit),
       y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.TempAnimalTypeSex[[m]] <- plots.TempAnimalTypeSex[[m]] +
    coord_cartesian(ylim = c(plots.TempAnimalTypeSex.limits$min[m], plots.TempAnimalTypeSex.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempAnimalTypeSex", ".svg")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempAnimalTypeSex[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempAnimalTypeSex", ".png")),
         width = 16, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.TempAnimalTypeSex[[m]], scale = 1)
  
  # results table
  tbl.plotResults.TempAnimalTypeSex <- rbind(tbl.plotResults.TempAnimalTypeSex, data.frame(rep(listOutcome[m], times = nrow(listResults.TempAnimalTypeSex.means[[1]])), listResults.TempAnimalTypeSex.means[[m]]$tempLevel, listResults.TempAnimalTypeSex.means[[m]]$animalType,  listResults.TempAnimalTypeSex.means[[m]]$sex, listResults.TempAnimalTypeSex.means[[m]]$emmean, listResults.TempAnimalTypeSex.means[[m]]$SE, listResults.TempAnimalTypeSex.means[[m]]$upper.CL - listResults.TempAnimalTypeSex.means[[m]]$emmean))
}

colnames(tbl.plotResults.TempAnimalTypeSex) <- c("targetVar", "tempLevel", "animalType", "sex", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TempAnimalTypeSex, file.path(pathWD, paste0("tblPlotResults", "_TempAnimalTypeSex", ".csv")), row.names = FALSE)

plots.TempAnimalTypeSex

```


# sex | animal type * temp, results
```{r}
listResults.SexAnimalTypeTemp <- purrr::map(listModel.final,
                                       ~ emmeans(.x,
                                                 pairwise ~ sex | animalType * tempLevel,
                                                 lmer.df = "satterthwaite",
                                                 adjust = NULL,
                                                 lmerTest.limit = 1e6)
)

listResults.SexAnimalTypeTemp.means <- purrr::map(listResults.SexAnimalTypeTemp,
                                             ~ .x$emmeans %>% as.data.frame()
)
listResults.SexAnimalTypeTemp.means <- purrr::map(listResults.SexAnimalTypeTemp.means,
                                                   ~ within(.x, {
                                                     lower.SEM <- emmean - SE
                                                     upper.SEM <- emmean + SE})
)

listResults.SexAnimalTypeTemp.contrasts <- purrr::map(listResults.SexAnimalTypeTemp,
                                                 ~ .x$contrasts %>% as.data.frame()
)

rm(listResults.SexAnimalTypeTemp)
```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listResults.SexAnimalTypeTemp.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listResults.SexAnimalTypeTemp.contrasts[[m]] <- listResults.SexAnimalTypeTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1) + 0.05, times = length(levels.AnimalType) * length(levels.Temp)),
           xEnd     = rep(c(2) - 0.05, times = length(levels.AnimalType) * length(levels.Temp)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots sex | temp * animal type
```{r fig.height=16/2.54, fig.width=14.93/2.54}
plots.SexTempAnimalType.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.5, 0.5),
                                             max = c(45, 10, 10, 0.3, 0.3, 2, 5.5))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.SexAnimalTypeTemp.contrasts[[m]] <- listResults.SexAnimalTypeTemp.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(1)) * (plots.SexTempAnimalType.limits$max[m] - plots.SexTempAnimalType.limits$min[m]) + plots.SexTempAnimalType.limits$min[m], times = length(levels.AnimalType) * length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.SexTempAnimalType.limits$max[m] - plots.SexTempAnimalType.limits$min[m]))
}

plots.SexTempAnimalType <- list()

for (m in 1:length(listModel.final)) {
  
  plots.SexTempAnimalType[[m]] <- ggplot(data = listResults.SexAnimalTypeTemp.means[[m]])
  
  # define layout parameters
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 12, angle = 90),
          strip.text.x = element_text(size = 12),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(tempLevel), rows = vars(animalType), labeller = labeller(animalType = labels.AnimalType, tempLevel = labels.Temp))
  
  # create plot
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    geom_segment(data = listResults.SexAnimalTypeTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.SexAnimalTypeTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.SexTempAnimalType.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.TempUnit,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = title.Sex,
       y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.SexTempAnimalType[[m]] <- plots.SexTempAnimalType[[m]] +
    coord_cartesian(ylim = c(plots.SexTempAnimalType.limits$min[m], plots.SexTempAnimalType.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexTempAnimalType", ".svg")),
         width = 14.93, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.SexTempAnimalType[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexTempAnimalType", ".png")),
         width = 14.93, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.SexTempAnimalType[[m]], scale = 1)
}

plots.SexTempAnimalType

```


## all plots sex | animal type * temp
```{r fig.height=16/2.54, fig.width=14.6/2.54}
plots.SexAnimalTypeTemp.limits <- data.frame(min = c(15, 0, 0, 0, 0, 0.0, 0.0),
                                             max = c(45, 10, 10, 0.3, 0.3, 1.75, 4.3))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listResults.SexAnimalTypeTemp.contrasts[[m]] <- listResults.SexAnimalTypeTemp.contrasts[[m]] %>%
    mutate(yStart = rep((0.89 + 0.08*c(0)) * (plots.SexAnimalTypeTemp.limits$max[m] - plots.SexAnimalTypeTemp.limits$min[m]) + plots.SexAnimalTypeTemp.limits$min[m], times = length(levels.AnimalType) * length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.SexAnimalTypeTemp.limits$max[m] - plots.SexAnimalTypeTemp.limits$min[m]))
}

plots.SexAnimalTypeTemp <- list()
tbl.plotResults.SexAnimalTypeTemp <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.SexAnimalTypeTemp[[m]] <- ggplot(data = listResults.SexAnimalTypeTemp.means[[m]])
  
  # define layout parameters
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 12, angle = 90),
          strip.text.x = element_text(size = 12),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(animalType), rows = vars(tempLevel), labeller = labeller(animalType = labels.AnimalType, tempLevel = labels.TempUnit))
  
  # create plot
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    geom_bar(mapping = aes(x = sex,
                           y = emmean,
                           fill = sex),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             fill = rep(colorPalette.AnimalTypeSex, times = 3),
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    geom_errorbar(mapping = aes(x = sex,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    scale_x_discrete(labels = labels.Sex)
  
  # add p-values
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    geom_segment(data = listResults.SexAnimalTypeTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.SexAnimalTypeTemp.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  # plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
  #   scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.SexAnimalTypeTemp.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.AnimalType,
       caption = paste("Model:", format(listFormula.full[[m]])),
       x = title.Sex,
       y = listLabelOutcome[m])
  
  # change coordinates = "zoom in"
  plots.SexAnimalTypeTemp[[m]] <- plots.SexAnimalTypeTemp[[m]] +
    coord_cartesian(ylim = c(plots.SexAnimalTypeTemp.limits$min[m], plots.SexAnimalTypeTemp.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexAnimalTypeTemp", ".svg")),
         width = 14.6, height = 16, units = "cm", dpi = 600, device = "svg",
         plot = plots.SexAnimalTypeTemp[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_SexAnimalTypeTemp", ".png")),
         width = 14.6, height = 16, units = "cm", dpi = 600, device = "png",
         plot = plots.SexAnimalTypeTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.SexAnimalTypeTemp <- rbind(tbl.plotResults.SexAnimalTypeTemp, data.frame(rep(listOutcome[m], times = nrow(listResults.SexAnimalTypeTemp.means[[1]])), listResults.SexAnimalTypeTemp.means[[m]]$sex, listResults.SexAnimalTypeTemp.means[[m]]$animalType, listResults.SexAnimalTypeTemp.means[[m]]$tempLevel,  listResults.SexAnimalTypeTemp.means[[m]]$emmean,  listResults.SexAnimalTypeTemp.means[[m]]$SE, listResults.SexAnimalTypeTemp.means[[m]]$upper.CL - listResults.SexAnimalTypeTemp.means[[m]]$emmean))
}

colnames(tbl.plotResults.SexAnimalTypeTemp) <- c("targetVar", "sex", "animalType", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.SexAnimalTypeTemp, file.path(pathWD, paste0("tblPlotResults", "_SexAnimalTypeTemp", ".csv")), row.names = FALSE)

plots.SexAnimalTypeTemp

```


# effect size, reference levels
```{r}
refLevel.effectSize.temp <- levels.Temp[1] # "Warm"

# labels for outcome variables: effect size
listLabelEffectSize <- c("Relative cell size",
                         expression("Relative " ~ tau["inlet"]),
                         expression("Relative " ~ tau["channel"]),
                         expression("Relative " ~ hat(d)["inlet"]),
                         expression("Relative " ~ hat(d)["channel"]),
                         "Relative Young's modulus",
                         "Relative viscosity")

```


# effect size: temp | animalType, results
```{r}
listEffectSizes.TempAnimalType.means <- list()
listEffectSizes.TempAnimalType.contrasts <- list()

for (m in 1:length(listModel.final)) {
  # copy contrast data
  listEffectSizes.TempAnimalType.contrasts[[m]] <- data.frame(contrast = listResults.TempAnimalType.contrasts[[m]]$contrast,
                                                              animalType = listResults.TempAnimalType.contrasts[[m]]$animalType,
                                                              df = listResults.TempAnimalType.contrasts[[m]]$df,
                                                              t.ratio = listResults.TempAnimalType.contrasts[[m]]$t.ratio,
                                                              p.value = listResults.TempAnimalType.contrasts[[m]]$p.value,
                                                              sig.idx = listResults.TempAnimalType.contrasts[[m]]$sig.idx,
                                                              p.valuef = listResults.TempAnimalType.contrasts[[m]]$p.valuef)
  
  # copy columns of tempLevel and animalType to results list for effect size
  listEffectSizes.TempAnimalType.means[[m]] <- data.frame(tempLevel = listResults.TempAnimalType.means[[m]]$tempLevel,
                                                          animalType = listResults.TempAnimalType.means[[m]]$animalType)
  # create columns for effect size and errors
  listEffectSizes.TempAnimalType.means[[m]]$mean <- NaN
  listEffectSizes.TempAnimalType.means[[m]]$SE <- NaN
  listEffectSizes.TempAnimalType.means[[m]]$lower.CL <- NaN
  listEffectSizes.TempAnimalType.means[[m]]$upper.CL <- NaN
  listEffectSizes.TempAnimalType.means[[m]]$lower.SEM <- NaN
  listEffectSizes.TempAnimalType.means[[m]]$upper.SEM <- NaN
  
  for (idx in 1:nrow(listEffectSizes.TempAnimalType.means[[m]])) {
    if (listEffectSizes.TempAnimalType.means[[m]]$tempLevel[idx] == refLevel.effectSize.temp) {
      # reference
      listEffectSizes.TempAnimalType.means[[m]]$mean[idx] <- 1
      listEffectSizes.TempAnimalType.means[[m]]$SE[idx] <- listResults.TempAnimalType.means[[m]]$SE[idx] / listResults.TempAnimalType.means[[m]]$emmean[idx]
      listEffectSizes.TempAnimalType.means[[m]]$lower.CL[idx] <- listResults.TempAnimalType.means[[m]]$lower.CL[idx] / listResults.TempAnimalType.means[[m]]$emmean[idx]
      listEffectSizes.TempAnimalType.means[[m]]$upper.CL[idx] <- listResults.TempAnimalType.means[[m]]$upper.CL[idx] / listResults.TempAnimalType.means[[m]]$emmean[idx]
      listEffectSizes.TempAnimalType.means[[m]]$lower.SEM[idx] <- listResults.TempAnimalType.means[[m]]$lower.SEM[idx] / listResults.TempAnimalType.means[[m]]$emmean[idx]
      listEffectSizes.TempAnimalType.means[[m]]$upper.SEM[idx] <- listResults.TempAnimalType.means[[m]]$upper.SEM[idx] / listResults.TempAnimalType.means[[m]]$emmean[idx]
    } else {
      # other level: calculate uncertainty by linear error term analysis
      # find reference data for same animal type
      refMask <- (listEffectSizes.TempAnimalType.means[[m]]$animalType == listEffectSizes.TempAnimalType.means[[m]]$animalType[idx]) & (listEffectSizes.TempAnimalType.means[[m]]$tempLevel == refLevel.effectSize.temp)
      # means of reference level and recent data row
      refLevel <- listResults.TempAnimalType.means[[m]]$emmean[refMask]
      dataLevel <- listResults.TempAnimalType.means[[m]]$emmean[idx]
      
      # error = confidence level
      refDelta <- listResults.TempAnimalType.means[[m]]$upper.CL[refMask] - listResults.TempAnimalType.means[[m]]$emmean[refMask]
      dataDelta <- listResults.TempAnimalType.means[[m]]$upper.CL[idx] - listResults.TempAnimalType.means[[m]]$emmean[idx]
      # effect size
      effectSize <- dataLevel / refLevel
      effectSizeDelta <- effectSize * sqrt( (dataDelta / dataLevel)^2 + (refDelta / refLevel)^2 )
      # assign results
      listEffectSizes.TempAnimalType.means[[m]]$mean[idx] <- effectSize
      listEffectSizes.TempAnimalType.means[[m]]$lower.CL[idx] <- effectSize - effectSizeDelta
      listEffectSizes.TempAnimalType.means[[m]]$upper.CL[idx] <- effectSize + effectSizeDelta
      
      # error = SEM
      refDelta <- listResults.TempAnimalType.means[[m]]$SE[refMask]
      dataDelta <- listResults.TempAnimalType.means[[m]]$SE[idx]
      # effect size
      effectSizeDelta <- effectSize * sqrt( (dataDelta / dataLevel)^2 + (refDelta / refLevel)^2 )
      listEffectSizes.TempAnimalType.means[[m]]$SE[idx] <- effectSizeDelta
      listEffectSizes.TempAnimalType.means[[m]]$lower.SEM[idx] <- effectSize - effectSizeDelta
      listEffectSizes.TempAnimalType.means[[m]]$upper.SEM[idx] <- effectSize + effectSizeDelta
      
    }
  }
}
rm(refMask, refLevel, dataLevel, refDelta, dataDelta, effectSize, effectSizeDelta, idx)
```


## effect size: fitted values
```{r}
listEffectSizes.fittedValues <- listResults.fittedValues

for (m in 1:length(listModel.final)) {
  listEffectSizes.fittedValues[[m]] <- listEffectSizes.fittedValues[[m]] %>%
    mutate(fitted.normalized = NA)
  
  for (idx in 1:nrow(listEffectSizes.fittedValues[[m]])) {
    # if (listEffectSizes.fittedValues[[m]]$tempLevel[idx] == refLevel.effectSize.temp) {
    #   # reference
    #   listEffectSizes.fittedValues[[m]]$fitted.normalized[idx] <- 1
    # } else {
      # other level:
      # find reference data for same animal type
      refMask <- (listResults.TempAnimalType.means[[m]]$animalType == listEffectSizes.fittedValues[[m]]$animalType[idx]) & (listResults.TempAnimalType.means[[m]]$tempLevel == refLevel.effectSize.temp)
      # means of reference level and recent data row
      refLevel <- listResults.TempAnimalType.means[[m]]$emmean[refMask]
      
      listEffectSizes.fittedValues[[m]]$fitted.normalized[idx] <- listEffectSizes.fittedValues[[m]]$fitted[idx] / refLevel
    # }
  }
}

```


## plot preparation
```{r}
# adapt contrasts data set for plotting
n = nrow(listEffectSizes.TempAnimalType.contrasts[[1]])
for (m in 1:length(listModel.final)) {
  listEffectSizes.TempAnimalType.contrasts[[m]] <- listEffectSizes.TempAnimalType.contrasts[[m]] %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.AnimalType)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.AnimalType)),
           xText  = 0.5*(xStart+xEnd))
}
```


## all plots effect size: temp | animalType v2
```{r fig.height=10/2.54, fig.width=16/2.54}
plots.effectSize.TempAnimalType.limits <- data.frame(min = c(.85, 0.5, 0.5, 0.5, 0.5, 0.75, 0.45),
                                                     max = c(1.15, 2.3, 2.3, 1.7, 1.7, 2.5, 3.8))
# prepare dataframe for p-values
for (m in 1:length(listModel.final)) {
  listEffectSizes.TempAnimalType.contrasts[[m]] <- listEffectSizes.TempAnimalType.contrasts[[m]] %>%
    mutate(yStart = rep((0.85 + 0.08*c(0,1,0)) * (plots.effectSize.TempAnimalType.limits$max[m] - plots.effectSize.TempAnimalType.limits$min[m]) + plots.effectSize.TempAnimalType.limits$min[m], times = length(levels.Temp)),
           yEnd = yStart,
           yText = yStart + 0.015*(plots.effectSize.TempAnimalType.limits$max[m] - plots.effectSize.TempAnimalType.limits$min[m]))
}

plots.effectSize.TempAnimalType <- list()
tbl.plotResults.effectSize.TempAnimalType <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.final)) {
  
  plots.effectSize.TempAnimalType[[m]] <- ggplot(data = listEffectSizes.TempAnimalType.means[[m]])
  
  # define layout parameters
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
    facet_grid(cols = vars(animalType), labeller = labeller(animalType = labels.AnimalType))
  
  # create point plot: replicate means
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    geom_beeswarm(data = listEffectSizes.fittedValues[[m]],
                  mapping = aes(x = tempLevel,
                                y = fitted.normalized),
                  colour = "#C0C0C0",
                  size = 2,
                  cex = 1.6)
  
  # create LMM plot: mean
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    geom_point(mapping = aes(x = tempLevel,
                             y = mean),
               stat = "identity", # important since we provide mean values already
               size = 3,
               colour = colorPalette.TempAnimalType)
  
  # add LMM plot: error bars
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  colour = colorPalette.TempAnimalType,
                  position = position_dodge())
  
  # add nested x-axis
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    geom_segment(data = listEffectSizes.TempAnimalType.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listEffectSizes.TempAnimalType.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, plots.effectSize.TempAnimalType.limits$max[m]),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.AnimalType,
         caption = paste("Model:", format(listFormula.full[[m]])),
         x = title.TempUnit,
         y = listLabelEffectSize[m])
  
  # change coordinates = "zoom in"
  plots.effectSize.TempAnimalType[[m]] <- plots.effectSize.TempAnimalType[[m]] +
    coord_cartesian(ylim = c(plots.effectSize.TempAnimalType.limits$min[m], plots.effectSize.TempAnimalType.limits$max[m]))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_effectSize", "_TempAnimalType", ".svg")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "svg",
         plot = plots.effectSize.TempAnimalType[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_effectSize", "_TempAnimalType", ".png")),
         width = 16, height = 10, units = "cm", dpi = 600, device = "png",
         plot = plots.effectSize.TempAnimalType[[m]], scale = 1)
  
  # results table
  tbl.plotResults.effectSize.TempAnimalType <- rbind(tbl.plotResults.effectSize.TempAnimalType, data.frame(rep(listOutcome[m], times = nrow(listEffectSizes.TempAnimalType.means[[1]])), listEffectSizes.TempAnimalType.means[[m]]$tempLevel, listEffectSizes.TempAnimalType.means[[m]]$animalType, listEffectSizes.TempAnimalType.means[[m]]$mean, listEffectSizes.TempAnimalType.means[[m]]$SE, listEffectSizes.TempAnimalType.means[[m]]$upper.CL - listEffectSizes.TempAnimalType.means[[m]]$mean))
}

colnames(tbl.plotResults.effectSize.TempAnimalType) <- c("targetVar", "tempLevel", "animalType", "mean", "SEM", "CI")
write.csv(tbl.plotResults.effectSize.TempAnimalType, file.path(pathWD, paste0("tblPlotResults", "_effectSize", "_TempAnimalType", ".csv")), row.names = FALSE)

plots.effectSize.TempAnimalType

```


# power dissipation: viscosity for temp (°C) | animal type
## calculations
```{r}
# A = 0.02939 mPa·s, B = 507.88 K, and C = 149.3 K (from Viswanath, D.S.; Natarajan, G. (1989). Data Book on the Viscosity of Liquids. Hemisphere Publishing Corporation. ISBN 0-89116-778-1, pp. 714–715.)
T <- c(10, 23, 37) + 273.15
A <- 0.02939 # mPa·s
B <- 507.88 # K
C <- 149.3 # K
viscWater_temp <- A * exp( B / (T - C) )
viscWater_temp
```


## plot preparations
```{r}
eta.hb.human <- 0.01 # Pa s
eta.hb.bat <- 0.1 * eta.hb.human # Pa s
index.eta <- which(listOutcome == "viscosity")
title.powerDiss <- "Power dissipation ratio\n(Cytoplasm / Membrane)"
colorPalette.powerDiss <- c("#C21618","#F3BA29","#1f77bd")

listResults.pd.TempAnimalType.means <- listResults.TempAnimalType.means[[index.eta]] %>%
  mutate(temp = case_when(tempLevel == "Cold" ~ 10,
                          tempLevel == "RT"   ~ 23,
                          tempLevel == "Warm" ~ 37,
                          TRUE ~ NaN)) %>%
  mutate(eta.hb = case_when(animalType == "Human" ~ eta.hb.human,
                            TRUE ~ eta.hb.bat)) %>%
  mutate(eta.hb = case_when(tempLevel == "Cold" ~ eta.hb * 1.3/.93,
                            tempLevel == "RT"   ~ eta.hb * .93/.93,
                            tempLevel == "Warm" ~ eta.hb * .69/.93,
                            TRUE ~ NaN)) %>%
  mutate(pd.ratio = eta.hb / emmean) %>%
  mutate(pd.ratio.upper.SEM = pd.ratio + eta.hb / emmean^2 * SE,
         pd.ratio.upper.SEM.2 = eta.hb / lower.SEM,
         pd.ratio.lower.SEM = pd.ratio - eta.hb / emmean^2 * SE,
         pd.ratio.lower.SEM.2 = eta.hb / upper.SEM) %>%
  mutate(temp = case_when(animalType == "RousettusAegyptiacus" ~ temp + .8,
                          animalType == "NyctalusNoctula" ~ temp - .8,
                          TRUE ~ temp))

```

## all plots power dissipation
```{r fig.height=10/2.54, fig.width=12/2.54}
plots.powerDiss.TempAnimalType.limits <- data.frame(min = c(1e-4),
                                                    max = c(2e-2))

tbl.plotResults.powerDiss.TempAnimalType <- data.frame(matrix(ncol = 5, nrow = 0))

plots.powerDiss.TempAnimalType <- ggplot(data = listResults.pd.TempAnimalType.means)

# exponential labels
fancy_scientific <- function(l) {
     # turn in to character string in scientific notation
     l <- format(l, scientific = TRUE)
     # quote the part before the exponent to keep all the digits
     l <- gsub("^(.*)e", "'\\1'e", l)
     # turn the 'e+' into plotmath format
     l <- gsub("e", "%*%10^", l)
     # return this as an expression
     parse(text=l)
}

# define layout parameters
plots.powerDiss.TempAnimalType <- plots.powerDiss.TempAnimalType +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.position = "none")

# create LMM plot: mean
plots.powerDiss.TempAnimalType <- plots.powerDiss.TempAnimalType +
  geom_point(mapping = aes(x = temp,
                           y = pd.ratio,
                           colour = animalType),
             stat = "identity", # important since we provide mean values already
             size = 3)

# add LMM plot: error bars
plots.powerDiss.TempAnimalType <- plots.powerDiss.TempAnimalType +
  geom_errorbar(mapping = aes(x = temp,
                              ymin = pd.ratio.lower.SEM,
                              ymax = pd.ratio.upper.SEM,
                              colour = animalType),
                size = .75,
                width = 1.2,
                position = position_dodge())

# add discrete x-axis labels
plots.powerDiss.TempAnimalType <- plots.powerDiss.TempAnimalType +
  scale_x_reverse(breaks = c(37,23,10))

# add further layout parameter
plots.powerDiss.TempAnimalType <- plots.powerDiss.TempAnimalType +
  scale_y_continuous(limits = c(plots.powerDiss.TempAnimalType.limits$min, plots.powerDiss.TempAnimalType.limits$max),
                     expand = expansion(mult = c(0.025, 0.025)),
                     trans = 'log10',
                     labels = fancy_scientific) +
  # scale_y_log10(limits = c(plots.powerDiss.TempAnimalType.limits$min, plots.powerDiss.TempAnimalType.limits$max),
  #               expand = expansion(mult = c(0.025, 0.025)),
  #               breaks = scales::trans_breaks("log10", function(x) 10^x),
  #               labels = scales::trans_format("log10", scales::math_format(10^.x))) +
  annotation_logticks(scaled = TRUE, sides = "l") +
  labs(caption = paste("Model:", format(listFormula.full)),
       x = paste(title.Temp, "(°C)"),
       y = title.powerDiss) +
  scale_color_manual(values = colorPalette.powerDiss,
                     breaks = names(labels.AnimalType.OneLine),
                     labels = labels.AnimalType.OneLine)

# change coordinates = "zoom in"
plots.powerDiss.TempAnimalType <- plots.powerDiss.TempAnimalType +
  coord_cartesian(ylim = c(plots.powerDiss.TempAnimalType.limits$min[m], plots.powerDiss.TempAnimalType.limits$max[m]))

# save plots
ggsave(filename = file.path(pathWD, paste0("powerDissipationRatio", "_TempAnimalType", ".svg")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "svg",
       plot = plots.powerDiss.TempAnimalType, scale = 1)

ggsave(filename = file.path(pathWD, paste0("powerDissipationRatio", "_TempAnimalType", ".png")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "png",
       plot = plots.powerDiss.TempAnimalType, scale = 1)

# # results table
tbl.plotResults.powerDiss.TempAnimalType <- rbind(tbl.plotResults.powerDiss.TempAnimalType, data.frame(rep("powerDissipationRatio", times = nrow(listResults.pd.TempAnimalType.means)), listResults.pd.TempAnimalType.means$tempLevel, listResults.pd.TempAnimalType.means$temp, listResults.pd.TempAnimalType.means$animalType, listResults.pd.TempAnimalType.means$pd.ratio, listResults.pd.TempAnimalType.means$pd.ratio.upper.SEM - listResults.pd.TempAnimalType.means$pd.ratio))


colnames(tbl.plotResults.powerDiss.TempAnimalType) <- c("targetVar", "tempLevel", "temp", "animalType", "mean", "SEM")
write.csv(tbl.plotResults.powerDiss.TempAnimalType, file.path(pathWD, paste0("tblPlotResults", "_powerDissipationRatio", "_TempAnimalType", ".csv")), row.names = FALSE)

plots.powerDiss.TempAnimalType

```


# Characteristic time for temp (°C) | animal type
```{r}
index.tau2 <- which(listOutcome == "tau2")
title.tau2 <- bquote(atop("Viscoelasticity","(Recovery time" ~ tau["c"] ~ "(ms))"))

listResults.tau2.TempAnimalType.means <- listResults.TempAnimalType.means[[index.tau2]] %>%
  mutate(temp = case_when(tempLevel == "Cold" ~ 10,
                          tempLevel == "RT"   ~ 23,
                          tempLevel == "Warm" ~ 37,
                          TRUE ~ NaN)) %>%
  mutate(temp = case_when(animalType == "RousettusAegyptiacus" ~ temp + 0.8,
                          animalType == "NyctalusNoctula" ~ temp - 0.8,
                          TRUE ~ temp))
```

## all plots tau2
```{r fig.height=10/2.54, fig.width=12/2.54}
plots.tau2.TempAnimalType.limits <- data.frame(min = c(1),
                                               max = c(4.5))

tbl.plotResults.tau2.TempAnimalType <- data.frame(matrix(ncol = 5, nrow = 0))

plots.tau2.TempAnimalType <- ggplot(data = listResults.tau2.TempAnimalType.means)

# define layout parameters
plots.tau2.TempAnimalType <- plots.tau2.TempAnimalType +
  theme_bw() +
  theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
        axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
        axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
        axis.text.y = element_text(size = 10, face = "bold"),
        panel.border = element_rect(size = 1.5, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
        plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
        plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
        legend.justification = c(0,1),
        legend.position = c(0.02,1.02),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.text.align = 0)

# create LMM plot: mean
plots.tau2.TempAnimalType <- plots.tau2.TempAnimalType +
  geom_point(mapping = aes(x = temp,
                           y = emmean,
                           colour = animalType),
             stat = "identity", # important since we provide mean values already
             size = 3)

# add LMM plot: error bars
plots.tau2.TempAnimalType <- plots.tau2.TempAnimalType +
  geom_errorbar(mapping = aes(x = temp,
                              ymin = lower.SEM,
                              ymax = upper.SEM,
                              colour = animalType),
                size = 0.75,
                width = 1.2,
                position = position_dodge())

# add discrete x-axis labels
plots.tau2.TempAnimalType <- plots.tau2.TempAnimalType +
  scale_x_reverse(breaks = c(10,23,37))

# add further layout parameter
plots.tau2.TempAnimalType <- plots.tau2.TempAnimalType +
  scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                     limits = c(0, plots.tau2.TempAnimalType.limits$max),
                     expand = expansion(mult = c(0.025, 0.025))) +
  labs(caption = paste("Model:", format(listFormula.full)),
       x = paste(title.Temp, "(°C)"),
       y = title.tau2) +
  scale_color_manual(values = colorPalette.powerDiss,
                     breaks = names(labels.AnimalType.OneLine),
                     labels = labels.AnimalType.OneLine)

# change coordinates = "zoom in"
plots.tau2.TempAnimalType <- plots.tau2.TempAnimalType +
  coord_cartesian(ylim = c(plots.tau2.TempAnimalType.limits$min, plots.tau2.TempAnimalType.limits$max))

# save plots
ggsave(filename = file.path(pathWD, paste0(listOutcome[index.tau2], "_v2", "_TempAnimalType", ".svg")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "svg",
       plot = plots.tau2.TempAnimalType, scale = 1)

ggsave(filename = file.path(pathWD, paste0(listOutcome[index.tau2], "_v2", "_TempAnimalType", ".png")),
       width = 12, height = 10, units = "cm", dpi = 600, device = "png",
       plot = plots.tau2.TempAnimalType, scale = 1)

# results table
tbl.plotResults.tau2.TempAnimalType <- rbind(tbl.plotResults.tau2.TempAnimalType, data.frame(rep(listOutcome[index.tau2], times = nrow(listResults.tau2.TempAnimalType.means)), listResults.tau2.TempAnimalType.means$tempLevel, listResults.tau2.TempAnimalType.means$temp, listResults.tau2.TempAnimalType.means$animalType, listResults.tau2.TempAnimalType.means$emmean, listResults.tau2.TempAnimalType.means$SE))

colnames(tbl.plotResults.tau2.TempAnimalType) <- c("targetVar", "tempLevel", "temp", "animalType", "mean", "SEM")
write.csv(tbl.plotResults.tau2.TempAnimalType, file.path(pathWD, paste0("tblPlotResults", "_v2", "_TempAnimalType", ".csv")), row.names = FALSE)

plots.tau2.TempAnimalType
```


# save data space image
```{r}
save.image(file.path(pathWD, paste0("RDataSpace_v8", ".RData")))
```

